#ifndef _SPARSE_MATRIX_H
#define _SPARSE_MATRIX_H

#include <iostream>
#include <strings.h>
#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>

class SparseMatrix
{
    friend std::ostream& operator<<(std::ostream&, const SparseMatrix&);

public:
    SparseMatrix (int rows, int columns, long active_cells)
        : m_rows(rows), m_active(active_cells), m_columns(columns), 
          m_curr_row(-1), m_curr_index(0)
    {
//        std::cout << "SparseMatrix(" << m_rows << "," << m_columns << "," 
//                  << m_active << ")" << std::endl;
        assert(m_columns > 0);
        assert(m_rows > 0);
        m_col_ind = new int [m_active];
        m_val = new double [m_active];
        m_row_ptr = new int [m_rows+1];
        m_row_ptr[m_rows] = m_active;
    }

    ~SparseMatrix ()
    {
        delete m_col_ind;
        delete m_row_ptr;
        delete m_val;
    }

    // matrix cells must be set in row/column order
    void set(int row, int column, double val)
    {
        assert(m_curr_row < m_rows);
        assert(column < m_columns);
        assert(row < m_rows);
        assert(row >= m_curr_row);

		for (int i=m_curr_row; i<row; ++i)
			m_row_ptr[++m_curr_row] = m_curr_index;

        m_col_ind[m_curr_index] = column; 
        m_val[m_curr_index] = val; 
        m_curr_index++;
        
        return;
    }
	// finishes off an unfinised row array
    void freeze()
	{
		for (int i=m_curr_row; i<m_rows; ++i)
			m_row_ptr[++m_curr_row] = m_curr_index;
	}

    boost::shared_array<double> 
	vector_product(boost::shared_array<double> vector, int start_row=0) const
    {
        boost::shared_array<double> result(new double [m_rows-start_row]);
        for(int i=start_row; i < m_rows; ++i) 
        {
            result[i-start_row] = 0;
            for(int j=m_row_ptr[i]; j < m_row_ptr[i+1]; ++j)
               result[i-start_row] += (m_val[j]*vector[m_col_ind[j]]);
        }
        return result;
    }

    double 
	vector_product_row(double *vector, int row) const
    {
//		std::cerr << row << std::endl;
		assert(row < m_rows);

        double result=0.0f;
		for(int j=m_row_ptr[row]; j < m_row_ptr[row+1]; ++j)
		   result += m_val[j]*vector[m_col_ind[j]];

        return result;
    }

    void
	get_row(long row, int &length, int **features, double **values)
    {
		assert(row < m_rows);
		*features=m_col_ind+m_row_ptr[row];
		*values=m_val+m_row_ptr[row];
		length=m_row_ptr[row+1]-m_row_ptr[row];;
    }

    boost::shared_array<double> 
	vector_product_transpose(boost::shared_array<double> vector, int start_row=0, 
							int end_row=0, bool mask=false) const
    {
        if (end_row == 0) end_row = m_rows;
        assert(start_row <= end_row);
        boost::shared_array<double> result(new double [m_columns]);
        bzero(result.get(),sizeof(double)*m_columns);

        if (mask)
        {
            for(int j=start_row; j < end_row; ++j) 
               for(int i=m_row_ptr[j]; i < m_row_ptr[j+1]; ++i)
                   result[m_col_ind[i]] += vector[j-start_row];
        }
        else
        {
            for(int j=start_row; j < end_row; ++j) 
            {
               for(int i=m_row_ptr[j]; i < m_row_ptr[j+1]; ++i)
                   result[m_col_ind[i]] += (m_val[i]*vector[j-start_row]);
            }
        }
        return result;
    }

    int rows( void ) const { return m_rows; }
    int columns( void ) const { return m_columns; }
    long active( void ) const { return m_active; }
    size_t size( void ) const 
    { return sizeof(double)*m_active + sizeof(int)*(m_active+m_rows); }

public:
    int *m_row_ptr;
    int *m_col_ind;
    double *m_val;
    
private:
    int m_rows; 
	long m_active;
    int m_columns;
    int m_curr_row, m_curr_index;
};
typedef boost::shared_ptr<SparseMatrix> SparseMatrixPtr;

inline std::ostream& operator<<(std::ostream &os, const SparseMatrix &sm)
{
    os << "<--------->" << std::endl;
    os << "Rows: " << sm.m_rows << " Columns: " << sm.m_columns 
		<< " Active: " << sm.m_active << std::endl;
    for(int row=0; row < sm.m_rows; ++row)
    {
        for(int it=sm.m_row_ptr[row]; it < sm.m_row_ptr[row+1]; ++it)
            os << sm.m_col_ind[it] << ":" << sm.m_val[it] << " ";
        os << std::endl;
    }
    os << "<--------->" << std::endl;

    return os;
}

#endif // _SPARSE_MATRIX_H
